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Abstract 

A  simple,  consistent,  and  transportable  magnitude  scale  for  regional  phases  is  useful  and  often 
required  in  order  to  (A)  improve  the  source  discrimination  capability,  (B)  determine  the  station  detection 
threshold,  and  (C)  estimate  the  explosive  yield.  A  convenient,  and  hence  recommended,  magnitude  for¬ 
mula  for  Lg  phase  is:^ 

.  4.0272  -  Bias  4  logA(A)  4  Itogwm))  4  l,og[sln(,j^^^f!^)l  4  ,  [1] 

where  the  “bias"  term  is  meant  to  account  for  the  different  Lg  excitation  relative  to  mt, .  For  instance,  a 
bias  of  approximately  0.39  magnitude  unit  for  Iranian  Plateau  has  been  suggested  by  Nuttli.® 

Given  a  suite  of  events  with  Lg  phases  recorded  at  a  seismic  network,  we  present  an  iterative  pro¬ 
cedure  to  simultaneously  invert  for  the  path  7  and  the  event  rrii^  values  in  Equation  [1]  without  using 
a  priori  path  7  information.  Independently  derived  7  (or  Q)  values,  if  available,  can  be  utilized  to  further 
constrain  the  trade-off  between  the  bias  term  in  [1]  and  the  resulting  7  values.  Other  constraints  can  be 
easily  incorporated  into  this  iterative  inversion  scheme  as  well.  The  procedure  is  less  sensitive  to  round¬ 
ing  errors,  and  hence  it  is  numerically  more  accurate  than  those  direct  methods  based  on  matrix  factori¬ 
zation  or  Gaussian  elimination.  When  the  number  of  equations  becomes  large,  the  iterative  approach  is 
often  the  only  practical  means  to  tackle  the  inversion.  This  joint  inversion  scheme  is  a  natural  extension 
to  crustal  phases  of  the  one  we  previously  used  in  the  teleseismic  analyses.** 

The  proposed  joint  inversion  scheme  has  been  tested  with  Pahute  Mesa  and  Novaya  Zemlya 
explosions,  and  the  -Lg  bias  at  these  two  sites  are  inferred  to  be  -0.34  and  -0.26  magnitude  unit, 
respectively.  These  bias  estimates  are  “optimal”  in  that  the  resulting  path  Oo  values,  which  are  the  by¬ 
product  of  this  inversion  exercise,  would  be  in  best  agreement  with  those  derived  by  the  coda-Q 
method.®  ®  ^  This  exercise  yields  twenty  and  eleven  stations  calibrated  for  Lg  phase  from  Pahute  Mesa 
and  Novaya  Zemlya  regions,  respectively. 

’  This  applied  research  was  solely  sponsored  by  the  Air  Force  Technical  Applications  Center.  Key  Words:  Lg,  magni¬ 
tude  scale,  bias,  transportability,  attenuation,  coda  Q,  inversion,  msiximum-likelihood  method. 
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Project  Objective 

Develop  simple,  transportable  magnitude  scales  for  miscellaneous  regional  phases.  Review  and 
improve  inversion  techniques  currently  in  use.  Some  anticipated  products/payoffs  include: 

[1]  a  more  consistent  way  of  calibrating  a  suite  of  propagation  paths, 

[2]  refined  event  magnitudes  which  have  immediate  or  potential  use  in  (a)  source  discrimination  study, 
(b)  determination  of  station  detection  threshold,  and  (c)  yield  estimation. 

Research  Accomplished 

To  date  the  absolute  Lg  magnitude  scale  is  defined  on  a  region-by-region  basis.  For  Eastern  U.S., 
Nuttli’s  (1986ab)  two-step  formulae  are  equivalent  to  the  following  one-step  procedure  (Jih  and  Lynnes, 
1993): 

m,.  .  4.0272  4  togA(A)  4  l|og(A(km))  4  4  ,  („ 

This  formula  defines  a  magnitude  scale  such  that  a  seismic  source  with  1-sec  Lg  amplitude  of  110 
pm  at  10  km  (extrapolated)  epicentral  distance  would  correspond  to  a  of  4.0272  +  2.0414  -i-  0.3333 
-  1.4019  +  0.0000  =  5.000,  which  was  suggested  to  be  appropriate  for  both  eastern  North  America  and 
Semipalatinsk.  That  is  to  say,  a  seismic  source  in  these  two  regions  with  5.0  would  have  a 
approximately  the  same. 

For  Iranian  Plateau,  Nuttll  (1980)  reported  that  seismic  sources  with  the  ISC  bulletin  5.0  excite 
Lg  amplitudes  approximately  270  microns  at  a  10-km  extrapolated  distance.  If  scale  is  to  be  “nor¬ 
malized”  to  rrit,  scale  at  =  5.0,  then  Equation  [1]  would  have  to  be  revised  for  Iran  as: 

.  4.0272  -  Bias  4  logA(A)  4  l|og(A(km))  4  l|og[sin(,pp;^fM^)l  4  .  [2] 

where  a  “bias”  of  approximately  0.39  magnitude  unit  [m.u.]  is  added  to  account  for  the  different  Lg  exci¬ 
tation  (relative  to  mt, )  observed  in  Iranian  Plateau.  Two  fundamental  issues  arise  immediately: 

[A]  If  we  accept  Equation  [2]  as  the  general  definition  of  iVi^  scale,  what  is  the  “bias”  term 
appropriate  for  other  places,  say  western  North  America  (such  as  NTS)  or  Novaya  Zemlya  regions? 

[B]  The  path  anelastic  attenuation  coefficient,  y,  is  assumed  to  be  known  before  Equation  [1]  (or 
[2])  can  be  applied.  How  should  the  path  attenuation  coefficient,  y,  be  determined  for  regions  like 
Novaya  Zemlya  where  Lg  might  be  blocked  rather  than  absorbed  through  the  intrinsic  attenuation? 
Would  the  coda-0  method  still  be  applicable  to  those  blocked  paths? 

Using  Lg  amplitudes  collected  under  a  recently  completed  AFTAC  contract  F08606-91-C-0005 
(Baumstark  and  Wagner,  1994),  these  two  issues  are  partially  examined  with  an  inversion  algorithm 
which  simultaneously  determines  the  event  sizes,  the  bias  term,  and  the  path  corrections  without  utiliz¬ 
ing  any  a  priori  path  y  information.  The  original  formulation  of  this  iterative  inversion  procedure  was  first 
presented  in  Jih  (1992),  and  it  was  tested  with  some  quasi-synthetic  data.  Its  updated  version  is  briefly 
described  in  the  Appendix  below. 

Given  a  postulated  bias  value,  a  system  of  linear  equations  (based  on  Equation  [2])  can  be  solved 
for  the  event  mi_^  and  path  y  values.  If  the  path  attenuation  coefficients  are  readily  available  from  other 
studies  independently,  then  such  extra  information  can  be  used  to  further  constrain  the  bias  term  for  the 


412 


17th  Seismic  Symposium 


11-15  Sept  1995 


most  probable  solution.  This  is  exactly  the  approach  to  be  used  in  the  following  exercises. 

l.  Novaya  Zemlya  Resuits 

For  Novaya  Zemlya  test  site,  there  are  24  explosions  recorded  at  13  stations,  totaling  92  Lg  paths. 
The  average  me  of  these  24  events  is  5.88,  based  on  the  path-corrected  /  station-corrected  me 
reported  in  Jih  and  Baumstark  (1994)  {cf.  pages  27-28).  This  value  was  used  to  constrain  the  inver¬ 
sion.  The  Qo  and  ii  values  associated  with  each  of  five  postulated  bias  terms  ranging  from  0.0  to  0.40 

m. u.  are  listed  in  Table  1.  It  turns  out  that  adopting  a  bias  of  0.26  m.u.  in  Equation  [2]  would  lead  to  the 
best  agreement  between  the  resulting  Qq  values  and  those  measured  by  Nuttli  (1988)  with  the  coda-Q 
method. 


Table  1 .  Oo ,  11  of  Novaya  Zemlya  -  WWSSN  Paths 

Station 

Nuttli 

Postulated  Lg  -m^  Bias  and  Resulting  Oo 

Code 

BSSA  1988 

.00 

.10 

.20 

.26* 

.40 

COP 

633  0.4 

802  0.86 

745  0.87 

697  0.89 

670  0.89 

615  0.91 

DAG 

290  0.68 

279  0.69 

268  0.71 

262  0.71 

249  0.73 

ESK 

499  0.67 

481  0.69 

463  0.70 

454  0.71 

433  0.73 

1ST 

569  0.72 

547  0.73 

536  0.74 

509  0.75 

KBS 

315  0.5 

KEV 

252  0.6 

314  0.50 

292  0.52 

274  0.54 

263  0.55 

242  0.57 

KON 

496  0.5 

454  0.04 

433  0.10 

421  0.13 

396  0.20 

243  0.47 

235  0.49 

228  0.51 

223  0.53 

213  0.55 

NUR 

420  0.5 

512  0.64 

479  0.66 

450  0.67 

434  0.68 

401  0.70 

STU 

531  0.5 

603  0.62 

577  0.64 

553  0.65 

540  0.66 

512  0.68 

TRI 

521  0.48 

502  0.51 

485  0.52 

475  0.54 

454  0.56 

UME 

391  0.5 

456  0.87 

427  0.88 

401  0.89 

386  0.89 

358  0.90 

2.  Pahute  Mesa  Results 

For  Pahute  Mesa  explosions,  a  Lg  -mt,  bias  of  0.34  m.u.  (Table  2)  appears  to  give  Qo  values 
most  consistent  with  those  Nuttli  (1986a)  and  Patton  (1988)  obtained.  225  Lg  signals  recorded  at  21 
Eurasian  stations  were  used  in  this  inversion.  The  47  Pahute  Mesa  events  have  an  average  ivg  of  5.51 
{cf.  pages  17-18  of  Jih  and  Baumstark,  1994). 
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Table  2.  C?o.  'H  of  Pahute  Mesa  -  WWSSN  Paths 

Station 

Nuttli+Patton 

Postulated  Lg  -mb  Bias  and  Resulting  Qq 

BSSA  86,  88 

.00 

.10 

.20 

.34* 

.40 

AAM 

904  1.11 

838  1.11 

782  1.10 

714  1.10 

689  1.10 

ALQ 

264  0.86 

247  0.87 

231  0.88 

212  0.89 

205  0.90 

ATL 

360  0.01 

352  0.07 

340  0.14 

336  0.17 

BKS 

139  0.6 

1800.88 

166  0.89 

153  0.89 

139  0.90 

133  0.90 

BLA 

560  0.38 

537  0.42 

515  0.45 

488  0.48 

476  0.50 

CMB 

138  0.36 

125  0.43 

115  0.49 

102  0.56 

98  0.59 

COR 

203  0.26 

193  0.30 

184  0.35 

173  0.39 

168  0.41 

ELK 

150  0.5 

287  0.47 

243  0.52 

211  0.56 

179  0.60 

168  0.61 

FVM 

373  0.12 

359  0.17 

346  0.21 

328  0.26 

321  0.28 

GOL 

234  0.53 

221  0.56 

209  0.58 

195  0.61 

190  0.62 

JAS 

149  0.66 

134  0.69 

123  0.71 

109  0.74 

104  0.75 

JCT 

456  0.72 

428  0.74 

402  0.76 

371  0.78 

359  0.79 

KNB 

142  0.4 

218  0.79 

181  0.81 

157  0.82 

132  0.82 

124  0.83 

LAC 

097  0.7 

189  0.53 

164  0.59 

145  0.65 

124  0.71 

117  0.73 

LON 

202  0.58 

194  0.60 

186  0.62 

176  0.65 

172  0.66 

LUB 

423  1.57 

393  1.55 

367  1 .52 

337  1 .49 

325  1 .48 

MNV 

093  0.6 

140  0.24 

110  0.45 

100  0.52 

OGD 

664  0.48 

637  0.52 

611  0.55 

577  0.58 

564  0.59 

520  0.04 

514  0.13 

506  0.20 

491  0.29 

484  0.32 

WES 

1205  1.02 

1114  1.01 

1035  1.01 

942  1 .00 

908  1.00 

Conclusions  and  Recommendations 

For  Pahute  Mesa,  where  the  Qq  (and  ti)  values  derived  by  the  coda-Q  method  are  believed  to  be 
appropriate  to  account  for  the  path  attenuation,  the  Lg-rrii,  bias  is  estimated  as  0.34  m.u.  This  value 
happens  to  be  in  agreement  with  the  published  bias  caused  by  the  different  upper  mantle  absorption 
in  the  eastern  and  western  U.S.  Thus  this  Pahute  Mesa  exercise  might  be  in  support  of  Nuttli’s  asser¬ 
tion  that  the  same  absolute  Lg  magnitude  scale  can  be  used  tor  both  eastern  and  western  U.S.  How¬ 
ever,  this  is  probably  an  exception  rather  than  a  general  rule. 

Previously  only  a  handful  number  of  stations  was  analyzed  by  Nuttli  (1986a)  and  Patton  (1988) 
with  the  coda-0  method.  We  now  have  20  “calibrated”  stations  for  Lg  phases  from  Pahute  Mesa.  How¬ 
ever,  the  proposed  calibration  is  subject  to  the  choice  of  the  “bias”  parameter  which  is  not  quite  obvious 
to  decide.  This  indeterminacy  of  Lg  absolute  magnitude  scale  could  be  a  persistent  issue  encountered 
in  every  attempt  of  using  seismic  phases  like  Lg  which  is  not  as  transportable  as  Ms  ■  A  (new)  feature 
of  this  simultaneous  inversion  code  is  that  it  offers  a  suite  of  solutions  to  choose  from.  If  only  a  “rela¬ 
tive”  Lg  scale  is  of  interest,  then  setting  the  bias  to  an  arbitrary  value,  say  0,  would  suffice.  In  any  case, 
the  simultaneous  inversion  can  calibrate  a  suite  of  stations  in  a  more  consistent  manner,  which  is  the 
typical  advantage  of  GLM  inversion  schemes. 

The  same  procedure  gives  a  Lg  -mb  bias  of  0.26  m.u.  for  northern  Novaya  Zemlya  explosions. 
Since  it  is  known  that  Lg  blockage  does  occur  for  paths  crossing  the  Barents  Shelf,  there  remains  a 
question  whether  it  is  appropriate  to  use  Nuttli’s  Qq  values  to  constrain  our  selection  of  the  bias  term. 


414 


17th  Seismic  Symposium 


11-15  Sept  1995 


Based  on  the  t'  study,  it  has  been  suggested  that  the  upper  mantle  of  Novaya  Zemlya  is  similar  to  that 
of  Semipalatinsk.  Thus  0.26  m.u.  is  expected  to  be  the  “upper  bound”  of  the  bias  term.  It  seems  that, 
however,  the  Oq  values  based  on  the  coda-0  method  can  provide  a  O  map  which  is  qualitatively  con¬ 
sistent  with  that  based  on  the  time-domain  computation.  The  Lg  blockage  at  Barents  Shelf  could  have 
caused  the  predominant  frequency  of  Lg  waves  to  shift  in  a  way  very  similar  to  that  due  to  a  stronger 
anelastic  attenuation.  Consequently,  Nuttli’s  Qg  values  may  be  biased  low,  yielding  an  over¬ 
compensation  of  the  attenuation  effect  in  computing  his  Novaya  Zemlya  iVi^  .  This  conjecture  can  be 
tested  with  numerical  modeling  experiments  using  LFD  method  {e.g.,  Jih,  1994). 

If  many  stations  are  deployed  along  the  same  azimuth  from  the  shot,  i.e. ,  if  a  profile  is  available, 
then  the  coda-Q  results  can  be  verified  independently  using  the  inter-station  spectral  ratios.  This  could 
be  included  in  future  field  experiments. 
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Appendix  A.  Joint  Inversion  Method:  Basic  Concepts 

Consider  the  case  of  m  explosions  recorded  at  some  or  all  of  n  stations.  For  simplicity,  we  assume 
that  these  m  explosions  are  detonated  in  the  same  test  site  for  the  moment.  The  case  of  multiple  test 
sites  will  be  discussed  later.  The  linear  model  for  Lg  phases  can  then  be  written  as 

E(i)  -  7(j)(A(i,j)-10km)/ln(10)  +  e(i,j)  =  Y(i.j)  ,  [3] 

where  Y(i,j)  =  4.0272  -  Bias  +  logA(i,j)  +  jlog(A(i,j))  +  ■^log[sin(  i^kivdegj'^^ ' 

Once  the  amplitudes  and  the  locations  of  the  events  (and  hence  the  epicentral  distances.  A)  are  avail¬ 
able,  Y  would  be  completely  known.  Only  the  event  sizes  (E)  and  the  path-specific  coefficients  of  ane- 
lastic  attenuation  (y)  are  the  unknown  parameters  to  be  determined.  The  obscuring  errors  e  are 
assumed  to  be  uncorrelated  and  to  belong  to  the  same  probability  distribution,  namely  a  common  Gaus¬ 
sian  distribution  with  zero  mean  and  variance  o^. 

If  all  the  events  are  clustered  in  a  small  region,  then  the  epicentral  distances  A(i,j)  would  be  almost 
identical  for  a  given  station  j.  In  this  case,  our  model  (Equation  [3])  is  a  special  case  of  a  more  general 
linear  system; 

E(i)  -(-  S(j)  +  e(i,j)  =  Y(i,j)  .  for  i  =  1 ,...,  m  ;  j  =  1 . n  .  [4] 

which  can  be  expressed  in  a  matrix  formulation: 

H  X  +  e  -  Y  ,  [5] 

where  H  is  the  design  (or  observation)  matrix.  X  and  Y  are  the  column  vectors  of  unknowns  and  obser¬ 
vations,  respectively.  The  standard  least-squares  [LS]  solution  (viz.  the  one  that  minimizes  the  residual 
sum  of  squares:  RSS  =  (Y  -  HX)'  (Y  -  HX))  to  any  linear  system  with  a  general  form  like  Equation  [5]  is 

Xls  s  (H'H)-^H'Y  .  [6] 

where  H'  is  the  transpose  of  H.  This  least-squares  estimator  has  many  optimality  properties.  For 
instance,  it  is  unbiased,  and  it  gives  minimum  variance  within  the  class  of  linear  unbiased  estimators. 
Furthermore,  ^ls  is  also  the  Maximum-Likelihood  Estimate  [MLE]  under  the  Gaussian  assumption.  It  is 
straightforward  to  compute  the  uncertainty  by  using  Var[XLs]  =  the  diagonal  of  a^(H'H)“\  which  is  simply 
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scaling  the  variance  of  the  random  perturbations  by  the  number  of  observations  associated  with  each 
unknown. 


In  our  case,  however,  the  matrix  H'H  in  Equation  [6]  is  singular,  and  hence  the  least-squares 
theory  can  not  be  applied  immediately  unless  the  linear  system  of  [4]  is  modified  somewhat.  Perhaps 
the  easiest  way  to  illustrate  the  indeterminacy  due  to  the  singularity  of  the  matrix  H'H  is  that  given  any 
set  of  solution  to  [4],  we  can  always  obtain  yet  another  set  of  solution  by  adding  a  constant  to  all  event 
magnitudes,  E(i),  i  =  1,...,  m,  and  subtracting  the  same  constant  from  each  station  term  S(j)  (Jih  and 
Shumway,  1989).  Alternatively,  we  can  verify  that  the  matrix  H'H  has  zero  determinant  with  linear  alge¬ 
bra  packages  such  as  UNPACK,  EISPACK,  or  LAPACK.  In  this  study,  however,  a  formal  proof  is 
presented  below.  Without  loss  of  generality,  we  can  assume  that  each  of  the  m  events  is  fully  recorded 
at  all  n  stations,  then 


X  =  [El,  Ea.  Ea,  .  .  .  ,  E^,  Si,  Sa, 


H'H  = 


n-ln 


"Imxn 

m-lp 


where  is  the  identity  matrix  of  order  m,  and  all  elements  of  the  m-by-n  matrix  l^xn  are  1.  For 
instance,  if  m  =  3  and  n  =  2,  then 


H'H  = 


2  0  0  1  1 

0  2  0  1  1 

0  0  2  1  1 

1113  0 
1110  3 


which  is  a  doubly-bordered  band  diagonal  sparse  matrix  (Tewarson,  1973;  Press  et  al.,  1988).  After 
exactly  n  row  operations  eliminating  the  lower-left  submatrix,  Inxm.  the  determinant  of  H'H  can  be  com¬ 
puted  (up  to  a  multiplicative  constant)  as  that  of 


Im  (•^)nxm 
Onxm  Pn 


m 


where  Pn  is  a  square  matrix  of  order  n  with  ^(n-1)  on  the  diagonal  and  —  elsewhere.  For  m  =  3  and 
n  =  2,  [7]  becomes 


1  0  0  0.5  0.5 

0  1  0  0.5  0.5 

0  0  1  0.5  0.5 

0  0  0  1.5  -1.5 

.0  0  0  -1.5  1.5 

It  suffices  to  examine  the  determinant  of  this  Pn,  or,  equivalently,  to  check  the  determinant  of  a  square 
matrix  of  order  n  with  n-1  on  the  diagonal  and  -1  elsewhere; 


n-1 

-1 

-1 

•  -1 

-1 

n-1 

-1 

•  -1 

-1 

-1 

n-1 

•  -1 

.-1 

-1 

-1 

•  n-1 

It  is  straightforward  to  prove  that,  by  mathematical  induction,  this  matrix  has  zero  determinant  for  any  n 
>  2.  Thus  the  matrix  H'H  in  our  linear  model  is  always  singular  regardless  how  good  the  observed 
amplitudes/magnitudes,  Y,  might  be.  We  therefore  need  an  extra  boundary  condition  to  constrain  our 
linear  model  for  a  unique  solution.  The  most  commonly  adopted  approach  in  teleseismic  magnitude 
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determination  is  to  have  all  station  terms  sum  up  to  zero.  This  implies  that  larger  events  would  tend  to 
be  unchanged  whether  we  apply  the  station  corrections  or  not,  and  hence  the  bulletin  magnitudes  of 
larger  events  published  by  ISC  and  NEIC  can  be  regarded  as  more  or  less  unbiased.  This  extra  con- 

n~1 

straint  can  be  incorporated  into  [4]  by  replacing  all  S{n)  by  -^SG)-  It  not  only  reduces  the  number  of 

unknowns  by  one,  but,  more  importantly,  regularizes  the  whole  linear  system  to  make  H'H  invertible. 
However,  this  is  not  the  only  plausible  constraint.  We  can  impose  the  extra  constraint  on  E(i)  instead,  or, 
even  impose  the  constraint  on  some  selected  stations. 

A  first  glance  of  Equation  [4]  might  lead  to  a  conclusion  that  the  inversion  scheme  for  crustal 
phases  is  identical  to  that  for  teleseismic  phase,  and  hence  the  algorithms  and  the  constraints  suitable 
for  the  teleseismic  data  reduction  would  be  appropriate  for  the  regional  case  as  well.  This  is  not  the 
case.  There  are  generally  fewer  regional  stations  available  for  inversions  with  crustal  phases.  The 
implicit  assumption  that  the  recording  stations  are  evenly  (and  randomly)  distributed  may  not  be  valid. 
Some  S(j)  terms  may  carry  more  weight  due  to  the  larger  corresponding  A(-,j).  This  implies  that  the 
zero-sum  assumption  on  the  S  terms  may  not  be  appropriate  for  the  i.^  inversion.  In  fact,  the  S  term  in 
Equation  [4],  by  definition,  can  not  have  zero  sum  across  the  network  because  both  y  and  A  in  Equation 
[3]  are  always  non-negative. 

Appendix  B.  Joint  Inversion  Method:  The  iterative  Procedure 

Once  a  constraint  has  been  chosen,  the  inverse  matrix  of  H'H  in  Equation  [6]  can  be  computed 
with  matrix  factorization  (such  as  Singular  Value  Decomposition,  [SVD])  or  Gaussian  elimination  method. 
Numerical  algorithms  of  these  types  are  called  direct  methods.  Direct  methods  can  be  impractical  if  H'H 
is  large.  In  that  case,  iterative  methods  are  often  the  only  possible  method  of  solution,  as  well  as  being 
faster  and  more  accurate  than  Gaussian  elimination  and  matrix  factorization.  The  largest  area  for  the 
application  of  iterative  methods  is  that  of  the  linear  systems  arising  in  the  numerical  solution  of  partial 
differential  equations.  Systems  of  orders  10,000  to  100,000  are  not  unusual  in  aerospace  sciences, 
although  the  majority  of  the  coefficients  of  the  systems  are  typically  zeros. 

The  basic  idea  of  iterative  methods  is  that  one  starts  with  a  trial  solution  vector  and  carries  out 
some  process  using  H,  Y,  and  X(°)  to  get  a  new  vector  X'^^  Then  one  repeats.  At  the  k  stage,  one  uses 
the  iterative  process  to  get  X^‘‘)  from  H  (or  H'H),  Y,  and  The  specific  algorithm  for  our  problem  is 
summarized  in  five  steps: 

Step  0 

Set  initial  value  of  yO)  for  j  =  1 . n. 

Step  1 

Compute  event  mi^  ,  E(i),  for  i  =  1 ,...,  m: 

E(i)  =  +  y(j)[A(i,j)-10km]/ln(10)] , 

where  #(])  is  the  number  of  stations  used  in  the  summation. 

Step  2 

Adjust  E(i),  i  =  1,...,  m,  with  the  desired  boundary  condition. 

Step  3 

Compute  the  path-specific  coefficient  of  anelastic  attenuation,  y(j),  for  j=1,...,  n: 
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7(j)  =  ln(10)B  E(i)  -  Y(i.j)  ]/BA{i,j)-10km] . 

i  i 


step  4 

Fill  in  missing  paths  (i,j)  with  predicted  pseudo-observations: 

Y(i.j)  =  E(i)-y(j)(A(i,j)-10km)/ln(10)  . 

Step  5 

Repeat  steps  [1]-[4]  to  update  E  and  y  till  convergence. 

Although  an  arbitrary  guess  of  y  will  do  at  Step  0,  picking  an  initial  y  close  to  the  average  across 
the  whole  area  of  interest  would  speed  up  the  convergence.  The  choice  of  the  extra  boundary  condition 
(at  Step  2)  is  very  flexible.  Constraining  the  mean  of  all  event  seems  to  perform  extraordinarily 
well,  however.  Note  that  the  pseudo-observations  predicted  at  Step  4  are  treated  as  “good”  observa¬ 
tions  at  steps  1  and  3  except  during  the  first  iteration  loop. 

The  iterative  procedure  described  above  is  a  special  case  of  a  general  algorithm  known  variously 
as  the  Gaussian-Seidel  method,  the  Liebmann  process,  or  the  method  of  successive  displacements 
(Bunch  and  Rose,  1976;  Forsythe  et  al.,  1977;  Golub  and  Van  Loan,  1983;  Spedicato,  1991).  The 
major  difference  between  this  method  and  that  of  Gaussian- Jacobi  (viz.,  the  simultaneous  displacement 
method)  is  that  we  solve  for  one  component  (viz,  E(i)  or  y{j))  of  the  new  vector  X  using  for  each  other 
component  of  X  its  most  recently  computed  value,  whereas  Gaussian-Jacobi  method  updates  all  unk¬ 
nown  parameters  simultaneously  at  the  end  of  each  iteration  loop.  In  our  case,  Gaussian-Seidel  method 
converges  faster  than  does  Gaussian-Jacobi  method.  The  first  application  of  this  iterative  technique  to 
the  magnitude  determination  problem  was  Blandford  and  Shumway  (1982),  although  the  methodology 
was  identified  as  the  E-M  [Expectation-Maximization]  algorithm  (Dempster  et  al.,  1977)  following  the 
convention  in  the  statistical  community. 

The  advantage  of  our  multi-event  joint  inversion  scheme,  as  compared  to  the  simple  network 
averaging  for  each  individual  event  that  Nuttli  (1986ab,  1988)  used,  is  that  we  can  have  a  more  con¬ 
sistent  network  for  all  events.  By  “consistent”  it  means  that  the  joint  inversion  procedure  provides  the 
best  approximation  to  the  network-averaged  values  that  would  have  been  obtained  if  all  the  events  were 
recorded  at  every  station  in  the  network.  This  advantage  may  not  be  very  obvious  in  using  direct 
methods  such  as  SVD  or  Gaussian  elimination.  However,  it  would  become  quite  natural  as  revealed  by 
Step  4  of  our  iteration  scheme. 
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